3. FEM-BEM Coupling#
The ngbem boundary element addon project initiated by Lucy Weggeler (see https://weggler.github.io/ngbem/intro.html) is now partly integrated into core NGSolve. Find a short and sweet introduction to the boundary element method there.
In this demo we simulate a plate capacitor on an unbounded domain.
from ngsolve import *
from netgen.occ import *
from ngsolve.solvers import GMRes, MinRes
from ngsolve.webgui import Draw
from ngsolve.bem import *
largebox = Box ((-2,-2,-2), (2,2,2) )
eltop = Box ( (-1,-1,0.5), (1,1,1) )
elbot = Box ( (-1,-1,-1), (1,1,-0.5))
largebox.faces.name = "outer" # coupling boundary
eltop.faces.name = "el1" # Dirichlet boundary
elbot.faces.name = "el2" # Dirichlet boundary
eltop.edges.hpref = 1
elbot.edges.hpref = 1
shell = largebox-eltop-elbot # FEM domain
shell.solids.name = "air"
mesh = shell.GenerateMesh(maxh=0.8)
mesh.RefineHP(2)
ea = { "euler_angles" : (-67, 0, 110) }
Draw (mesh, clipping={"x":1, "y":0, "z":0, "dist" : 1.1}, **ea);
On the exterior domain \(\Omega^c\), the solution can be expressed by the representation formula:
\[ x \in \Omega^c: \quad u(x) = - \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, \gamma_1 u (y)\, \mathrm{d}\sigma_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}\sigma_y\,, \]
where \(\gamma_0 u = u\) and \(\gamma_1 u = \frac{\partial u}{\partial n}\) are Dirichlet and Neumann traces. These traces are related by the Calderon projector
\[ \begin{align}\begin{aligned}\begin{split} \left( \begin{array}{c} \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{cc} -V & \frac12 + K \\ \frac12 - K^\intercal & -D \end{array} \right) \left( \begin{array}{c} \gamma_1 u \\ \gamma_0 u \end{array}\right) $$.\end{split}\\The $V$, $K$ are the single layer and double layer potential operators, and $D$ is the hypersingular operator.\\On the FEM domain we have the variational formulation\end{aligned}\end{align} \]
\int_{\Omega_\text{FEM}} \nabla u \nabla v , dx - \int_\Gamma \gamma_1 u v , ds = 0 \qquad \forall , v \in H^1(\Omega_\text{FEM}) $$
We use Calderon’s represenataion formula for the Neumann trace:
\[
\int_{\Omega_\text{FEM}} \nabla u \nabla v \, dx - \int_\Gamma \left( \left( \tfrac{1}{2} - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) v = 0 \qquad \forall \, v \in H^1(\Omega_\text{FEM})
\]
To get a closed system, we use also the first equation of the Calderon equations. To see the structure of the discretized system, the dofs are split into degrees of freedom inside \(\Omega\), and those on the boundary \(\Gamma\). The FEM matrix \(A\) is split accordingly. We see, the coupled system is symmetric, but indefinite:
\[\begin{split} \left( \begin{array}{ccc } A_{\Omega\Omega} & A_{\Omega\Gamma} & 0 \\ A_{\Gamma\Omega} & A_{\Gamma\Gamma } + D & -\frac12 M^\intercal + K^\intercal \\ 0 & -\frac12 M + K & -V \end{array}\right) \left( \begin{array}{c} u \\ \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{c} F_{\Omega} \\ F_{\Gamma}\\ 0 \end{array}\right) \,. \end{split}\]
Generate the finite element space for \(H^1(\Omega)\) and set the given Dirichlet boundary conditions on the surfaces of the plates:
order = 4
fesH1 = H1(mesh, order=order, dirichlet="el1|el2")
print ("H1-ndof = ", fesH1.ndof)
H1-ndof = 91269
The finite element space \(\verb-fesH1-\) provides \(H^{\frac12}(\Gamma)\) conforming element to discretize the Dirichlet trace on the coupling boundary \(\Gamma\). However we still need \(H^{-\frac12}(\Gamma)\) conforming elements to discretize the Neumann trace of \(u\) on the coupling boundary. Here it is:
fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries("outer"))
print ("L2-ndof = ", fesL2.ndof)
L2-ndof = 4035
fes = fesH1 * fesL2
u,dudn = fes.TrialFunction()
v,dvdn = fes.TestFunction()
a = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble()
gfudir = GridFunction(fes)
gfudir.components[0].Set ( mesh.BoundaryCF( { "el1" : 1, "el2" : -1 }), BND)
f = LinearForm(fes).Assemble()
res = (f.vec - a.mat * gfudir.vec).Evaluate()
Generate the the single layer potential \(V\), double layer potential \(K\) and hypersingular operator \(D\):
n = specialcf.normal(3)
with TaskManager():
V = LaplaceSL(dudn*ds("outer"))*dvdn*ds("outer")
K = LaplaceDL(u*ds("outer"))*dvdn*ds("outer")
D = LaplaceSL(Cross(grad(u).Trace(),n)*ds("outer"))*Cross(grad(v).Trace(),n)*ds("outer")
M = BilinearForm(u*dvdn*ds("outer"), check_unused=False).Assemble()
Setup the coupled system matrix and the right hand side:
sym = a.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat
rhs = res
bfpre = BilinearForm(grad(u)*grad(v)*dx+1e-10*u*v*dx + dudn*dvdn*ds("outer") ).Assemble()
pre = bfpre.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
Compute the solution of the coupled system:
with TaskManager():
sol_sym = MinRes(mat=sym, rhs=rhs, pre=pre, tol=1e-8, maxsteps=200, printrates=True)
LinearSolver iteration 1, residual = 18.46865730865644
LinearSolver iteration 2, residual = 2.3400628305255675
LinearSolver iteration 3, residual = 2.082508582473063
LinearSolver iteration 4, residual = 0.42147022859368977
LinearSolver iteration 5, residual = 0.26625863647556525
LinearSolver iteration 6, residual = 0.112055192022602
LinearSolver iteration 7, residual = 0.05786199287911832
LinearSolver iteration 8, residual = 0.035293799810918056
LinearSolver iteration 9, residual = 0.015622853897268196
LinearSolver iteration 10, residual = 0.014077891549913601
LinearSolver iteration 11, residual = 0.01405343484261758
LinearSolver iteration 12, residual = 0.0035502477934946475
LinearSolver iteration 13, residual = 0.0030969662181524115
LinearSolver iteration 14, residual = 0.00208372462376382
LinearSolver iteration 15, residual = 0.0016137942053708325
LinearSolver iteration 16, residual = 0.0013415743220771112
LinearSolver iteration 17, residual = 0.001327577684814304
LinearSolver iteration 18, residual = 0.0007480874466254629
LinearSolver iteration 19, residual = 0.0007138749392613201
LinearSolver iteration 20, residual = 0.0005701981932483101
LinearSolver iteration 21, residual = 0.0005637140390841748
LinearSolver iteration 22, residual = 0.0003838635724185585
LinearSolver iteration 23, residual = 0.000383552206745772
LinearSolver iteration 24, residual = 0.0002819116996648131
LinearSolver iteration 25, residual = 0.0002735292097634044
LinearSolver iteration 26, residual = 0.00021513665389341754
LinearSolver iteration 27, residual = 0.00020073895252908615
LinearSolver iteration 28, residual = 0.00013795110719783604
LinearSolver iteration 29, residual = 0.00013767309266065432
LinearSolver iteration 30, residual = 0.00012064207932672481
LinearSolver iteration 31, residual = 0.00010761004074587245
LinearSolver iteration 32, residual = 7.815674568311012e-05
LinearSolver iteration 33, residual = 7.766309208130407e-05
LinearSolver iteration 34, residual = 6.965065129226461e-05
LinearSolver iteration 35, residual = 6.197480566073734e-05
LinearSolver iteration 36, residual = 6.110385489291538e-05
LinearSolver iteration 37, residual = 4.408912682452194e-05
LinearSolver iteration 38, residual = 4.003428991466895e-05
LinearSolver iteration 39, residual = 3.7069138281061434e-05
LinearSolver iteration 40, residual = 3.060671627932603e-05
LinearSolver iteration 41, residual = 2.726232413327262e-05
LinearSolver iteration 42, residual = 2.3118721033224532e-05
LinearSolver iteration 43, residual = 2.218783704796981e-05
LinearSolver iteration 44, residual = 2.106234376657882e-05
LinearSolver iteration 45, residual = 1.911298410467399e-05
LinearSolver iteration 46, residual = 1.907632343957952e-05
LinearSolver iteration 47, residual = 1.7497820405561993e-05
LinearSolver iteration 48, residual = 1.574882445805853e-05
LinearSolver iteration 49, residual = 1.2285897257211445e-05
LinearSolver iteration 50, residual = 1.223938339838736e-05
LinearSolver iteration 51, residual = 1.0054615220491016e-05
LinearSolver iteration 52, residual = 1.0036144027854566e-05
LinearSolver iteration 53, residual = 7.952985414567752e-06
LinearSolver iteration 54, residual = 7.952710283421163e-06
LinearSolver iteration 55, residual = 6.6621413165794125e-06
LinearSolver iteration 56, residual = 6.654792823255396e-06
LinearSolver iteration 57, residual = 6.503704369769871e-06
LinearSolver iteration 58, residual = 5.592331723154123e-06
LinearSolver iteration 59, residual = 4.776657520114305e-06
LinearSolver iteration 60, residual = 4.67624306879596e-06
LinearSolver iteration 61, residual = 4.376979340880257e-06
LinearSolver iteration 62, residual = 3.955551462840321e-06
LinearSolver iteration 63, residual = 3.6361790132357767e-06
LinearSolver iteration 64, residual = 3.248804953802332e-06
LinearSolver iteration 65, residual = 2.7792991013866356e-06
LinearSolver iteration 66, residual = 2.749410133923865e-06
LinearSolver iteration 67, residual = 2.4594865872596695e-06
LinearSolver iteration 68, residual = 2.369681628115071e-06
LinearSolver iteration 69, residual = 2.3565749408041453e-06
LinearSolver iteration 70, residual = 1.9582827922013903e-06
LinearSolver iteration 71, residual = 1.939002338425858e-06
LinearSolver iteration 72, residual = 1.589831938180481e-06
LinearSolver iteration 73, residual = 1.582106754153215e-06
LinearSolver iteration 74, residual = 1.4424918548409193e-06
LinearSolver iteration 75, residual = 1.4284464030284612e-06
LinearSolver iteration 76, residual = 1.0656347382967589e-06
LinearSolver iteration 77, residual = 1.0463456362932312e-06
LinearSolver iteration 78, residual = 9.060701177013008e-07
LinearSolver iteration 79, residual = 9.016499256309143e-07
LinearSolver iteration 80, residual = 8.350770755807928e-07
LinearSolver iteration 81, residual = 8.055004097388477e-07
LinearSolver iteration 82, residual = 6.266897100478454e-07
LinearSolver iteration 83, residual = 6.250647210715121e-07
LinearSolver iteration 84, residual = 5.001546337658163e-07
LinearSolver iteration 85, residual = 4.972518927319457e-07
LinearSolver iteration 86, residual = 4.178706010970862e-07
LinearSolver iteration 87, residual = 4.08775595385533e-07
LinearSolver iteration 88, residual = 3.372181202794672e-07
LinearSolver iteration 89, residual = 3.359103038183894e-07
LinearSolver iteration 90, residual = 3.134458366039041e-07
LinearSolver iteration 91, residual = 2.9239466773053397e-07
LinearSolver iteration 92, residual = 2.8955747721339424e-07
LinearSolver iteration 93, residual = 2.517829983043718e-07
LinearSolver iteration 94, residual = 2.46799835101669e-07
LinearSolver iteration 95, residual = 2.1464248290005083e-07
LinearSolver iteration 96, residual = 2.1039934311767813e-07
LinearSolver iteration 97, residual = 1.6642984779893303e-07
gfu = GridFunction(fes)
gfu.vec[:] = sol_sym + gfudir.vec
Draw(gfu.components[0], clipping={"x" : 1, "y":0, "z":0, "dist":0.0, "function" : True }, **ea, order=2);
The Neumann data:
Draw (gfu.components[1], **ea, order=3);